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Abstract 

We explore a framework to model the dose response of allosteric multisite phosphorylation proteins using a single auxiliary 
variable. This reduction can closely replicate the steady state behavior of detailed multisite systems such as the Monod- 
Wyman-Changeux allosteric model or rule-based models. Optimal ultrasensitivity is obtained when the activation of an 
allosteric protein by its individual sites is concerted and redundant. The reduction makes this framework useful for modeling 
and analyzing biochemical systems in practical applications, where several multisite proteins may interact simultaneously. 
As an application we analyze a newly discovered checkpoint signaling pathway in budding yeast, which has been proposed 
to measure cell growth by monitoring signals generated at sites of plasma membrane growth. We show that the known 
components of this pathway can form a robust hysteretic switch. In particular, this system incorporates a signal proportional 
to bud growth or size, a mechanism to read the signal, and an all-or-none response triggered only when the signal reaches a 
threshold indicating that sufficient growth has occurred. 
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Introduction 

Protein phosphorylation is a common form of post-translational 
modification frequently used in nature to alter protein activity, for 
instance by changing the electrostatic properties of the protein or 
its spatial structure. The phosphorylation of the same protein at 
multiple different aminoacid residues is also very common, and it 
is found in proteins such as p53 [1], Sicl [2], EGFR [3], Weel [4], 
Ste5 [5] and many others [6] . 

The differences in the function of single-site vs. multisite 
phosphorylation are not completely understood. Many multisite 
proteins are involved in regulatory processes that can benefit from 
the presence of bistability, hysteresis, or limit cycles, which require 
sufficiently nonlinear interactions in addition to the right type of 
feedback [7,8]. A reasonable hypothesis is that multisite phos- 
phorylation can give rise to ultrasensitive dose responses, in a way 
that would not be possible in a comparable single-site system 
[9-13]. Many detailed mechanisms have also been proposed to 
explain the role of multisite systems in the emergence of bistability 
(for examples, see [14,15]). 

On the other hand, such detailed multisite mechanisms are 
normally not used as part of actual mathematical models of 
biochemical interactions. This is because explicitly modeling 
multiple sites usually involves the introduction of numerous 
variables, one or more for each phosphorylation state, and realistic 
models are often too complex already to justify this additional effort. 
Systems that attempt to model biochemical reactions explicidy often 
use the assumption that the protein has two states, one active and 



one inactive, with a simple reaction to transform one into the other, 
effectively assuming that the protein only has one site. Other models 
are more phenomenological in nature and include, for example, 
Hill function terms in the equation that are less clearly tied to the 
actual biochemical reactions [16,17]. 

In this paper we describe a simple mechanistic approach for 
modeling multisite allosteric proteins. This approach, named 
modified fraction (MF) modeling, is capable of describing 
ultrasensitive dynamics without introducing a large number of 
additional variables. Under this framework one keeps track of the 
fraction of modified sites in the protein, and the concentration of 
active protein over time is estimated from this information. The 
protein is activated in a way that requires the phosphorylation of 
several but not all of the sites. The approximation becomes 
increasingly precise as the number of sites increases, with good 
estimates already for around four or more sites. In a sense, this 
mechanism can be considered a one-variable, quasi-steady state 
reduction of a model similar to the Monod-Wyman-Changeux 
allosteric system [18], although uses of MF modeling outside of 
MWC are also possible. The MF framework can also be extended 
to other types of multisite modification such as ligand binding, 
multisite transcription factor regulation, multisite methylation or 
acetylation, ubiquitination, etc [19,20]. 

Perhaps the best way to test the versatility of a modeling tool 
such as the one proposed is to implement it in an actual 
biochemical system. In the current work we describe a detailed 
mechanistic model of a cell size checkpoint. Cell size checkpoints 
halt the cell cycle at specific points until sufficient cell growth has 
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Author Summary 

A large number of proteins in the cell are modified post- 
translationally by phosphorylation at multiple specific 
locations. This can help bring about interesting dynamical 
behaviors such as bistability or all-or-none responses to 
stimuli. Such behaviors are in turn important for cellular 
decision-making, differentiation, or the regulation of 
cellular processes. In this paper we propose the use of a 
specific technique for modeling allosteric multisite pro- 
teins, which can be thought of as the reduced version of a 
more detailed set of reactions. After validating this 
technique computationally by comparison to more de- 
tailed systems, we apply it to a new model of a signal 
transduction pathway with several multisite proteins. The 
model is concerned with a mechanism in budding yeast 
that is thought to measure the extent of daughter bud 
growth and send a signal to initiate mitosis when sufficient 
growth has occurred. Using the given framework we 
derive an analytically tractable model that creates the 
desired all-or-none signal. Overall, we give quantitative 
support to a newly discovered biochemical pathway, and 
we show the utility of the new modeling framework in the 
context of a realistic biological problem. 



occurred [21,22]. The mechanisms by which cell size checkpoints 
operate are poorly understood, and it is unclear whether they 
monitor actual cell size or parameters more closely related to the 
extent or rate of growth. In budding yeast, growth of a new cell is 
initiated when a daughter bud is formed on the surface of the cell 
[23]. The daughter bud initially grows in a polar manner, with all 
growth directed to the bud tip. Growth of the bud eventually 
switches to isotropic growth, in which the bud grows over its entire 
surface (see Figure 1A) [24]. The timing of the switch determines 
the duration of polar growth, which influences cell size and shape. 
It has been proposed in recent work by one of the authors that a 
cell size checkpoint controls the timing of this switch [25]. This 
checkpoint is the subject of our model. The variables are 
illustrated in Figure IB and described in more detail below. See 
Tyson and Novak [26] for an accessible introduction to the 
systems-wide modeling of cell cycle checkpoints. 

The MF approximation equation was first developed in [27], in 
the context of multisite systems with independent modification 
sites, with an emphasis on the estimation of the Hill exponents of 
sequential and nonsequential systems and on the comparison of 
their qualitative behavior. A major advance of the current paper, 
beyond the application to the cell cycle checkpoint, is to extend 
this work to cooperative and allosteric systems. Such systems are 
by definition non-independent, since the modification of one site 
accelerates the rate of modification of its neighbors. Cooperative 
systems are also more common and much better characterized 
than independent ones. The validation of the approximation in 
cooperative systems is ultimately based on a computational 
comparison of the MF reduction with detailed cooperative models 
having n or 2" variables. 

In the first two Results sections we carry out a description of the 
modified fraction method to model multisite systems, and we 
compare simulations of the reduced model with those of a detailed 
mechanistic model. In the remaining two Results sections we carry 
out a mathematical analysis of the proposed checkpoint signaling 
pathway. We hypothesize that this interaction pathway has the 
capacity to produce a bistable signal responsible for a sudden 
switch from polar to isotropic growth, once the bud has undergone 
sufficient polar growth. The model presents several desirable 
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Figure 1. The Rhol Network. A: A yeast bud grows first in a 
particular direction (polar growth) and eventually switches to growth in 
all directions (isotropic growth). B: The Rho1 signaling pathway starts 
with inactive Rhol flowing into the bud attached to membrane vesicles. 
Rho1 is then activated and binds to Pkd, forming an upstream system. 
The downstream system describes the activation of the PP2A/Zds1 
dimer leading up to the modification of cell cycle regulatory protein 
Cdkl. Multiple intermediate feedback loops are shown to allow for 
robust hysteretic and switch-like behavior. 
doi:10.l371/journal.pcbi.l003443.g00l 

qualities for a checkpoint, in particular a clear downstream signal 
when a sufficient polar bud growth has occurred. 

Results 

The MF framework describes a compact multisite 
mechanism 

We start by describing the assumptions on our model in the 
context of multisite phosphorylation (although it could also be 
applied to other irreversible covalent modifications as well as 
noncovalent ligand binding). Suppose that a protein substrate P is 
phosphorylated by a kinase E at n possible sites and dephosphor- 
ylated by a phosphatase F. The system is assumed to be 
nonsequential, so that there are 2" different phosphoforms of P, 
and the number of sites is thought to be relatively large e.g. n>5. 
The system is cooperative in the sense that site phosphorylation 
accelerates the phosphorylation of neighboring sites. Since the 
number of sites is relatively large, the activation is thought to be 
cumulative and the effect of any individual site is assumed to be 
small. The sites are assumed to be equivalent to each other, in the 
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sense that the rate of phosphorylation and dephosphorylation is 
similar across all sites and that no site has a stronger effect on 
substrate activation than other sites. The activation of the substrate 
may be due to binding to another molecule or body, such as the 
cell membrane. It could also be due to an internal structural 
change that allows the substrate to interact differently with other 
proteins. Thus the active protein concentration can be denned as 
the concentration of the protein bound to a particular molecule or 
in a particular molecular state. 

Suppose that the phosphorylations lead to the concerted, 
redundant activation of the protein. That is, multiple phosphor- 
ylations are necessary for activation (concerted), and not all sites 
need to be phosphorylated for full activation (redundant). We 
define the activity function h : [0,1] — > [0,1] such that the fraction of 
active protein with i phosphorylated sites is given by h(i/n). 

There are many systems that likely fall within this general 
framework. For instance, Ste5 is a scaffold protein in budding 
yeast with n = 8 phosphorylation sites, which relays a pheromone 
response only when it is bound to the membrane [28,29] . When 
phosphorylated by Cdkl, it tends to unbind from the membrane, 
shutting down its activity. The sites are predicted to lie on an 
unstructured region of the protein and appear to act by changing 
the protein's bulk electrostatic properties. In the paper [29], it was 
shown through site-directed mutagenesis that around five or more 
phosphorylations are necessary and sufficient for deactivation. 
Another recent example is the multisite phosphorylation of Cdc25 
by Cdkl in fission yeast, which was similarly studied in detail by 
mutating individual sites [30]. 

According to the modified fraction framework, we estimate the 
concentration of a particular protein state from the overall fraction 
of modified sites. For instance, if the protein has w = 3 sites and the 
fraction of phosphorylated sites is p, then the fraction of protein 
with only the first and last sites phosphorylated is roughly 

101 2/i % 

Here P is the total protein concentration. This is not an equality 
since cooperative effects introduce correlations among the sites, i.e. 
the sites are not independent of each other, but it is an 
approximation assuming cooperative effects are sufficiently weak. 
Multiplying on both sides by P and adding over all possible 
phosphoforms with i phosphorylations, the concentration P, of 
proteins with exactly i phosphorylated sites out of a total of n sites is 

p^P^Ai-pf-. 

One can estimate the overall concentration P of active protein as 

A key aspect of this formula is that if we denote the right hand side 
by f n (p), it can be shown that f n (p) converges to Ph(p) as n increases 
[3 1] . This gives the approximation 

P~Ph{p), (1) 

which becomes increasingly precise for large n. Notice that the 
different quantities in this formula can potentially be measured in 
the lab - the active protein concentration via an activity assay, the 



total protein concentration via Western blot, and the activity 
function h(x) through site-directed mutagenesis. 

A timescale decomposition argument can be made to use this 
approximation away from steady state. If p(i) is the fraction of 
active sites over time, and the timescale of protein activation is 
much faster than the rate at which p(t) changes, then one can 
approximate P{t)xPh(p{t)) at any given time using the same 
formula. This produces a convenient method for modeling 
multisite systems under the given assumptions by keeping track 
of the variable p(j), without creating n, let alone 2" variables. On 
the other hand, if p(t), or any other process affecting protein 
activation, is at least as fast as protein activation itself, then 
nontrivial dynamics might take place such as limit cycle 
oscillations, and the approximation can introduce errors. 

It is necessary to calculate the fraction of phosphorylation p 
itself. Assuming linear rates of phosphorylation and dephosphor- 
ylation for p, one obtains the system 

«£ 

Pin T± P- 

fSF 

That is, dp I dt = nEp in — fSFp, where p in = 1 — p is the fraction of 
inactive sites. This is the default for the model, although any other 
rate equation for dp/ dt can be used, including Michaelis-Menten 
complex formation at the level of the individual sites. 

As for the activity function h(x), any sigmoidal function can be 
used, including functions measured directly by experiments. By 
default we assume the following form, which we will derive in the 
next section: 



See also [5,9,27] for other uses and derivations of this formula in 
the literature. The function h(x) can actually be highly switch-like 
for large y, which illustrates how small increases in the kinase E 
can result in large activity changes in the protein P, unlike linear 
rate models with only one site. Other forms for the activity 
function h(x) have effectively been considered by by Kapuy et al 
[32], and also by Wang et al [10]. 

In Figure 2A we show the relationship between the fraction of 
phosphorylated sites p and the active protein concentration P for a 
particular choice of the MF parameters tX,P,y,d, using the 
approximation formula PxPh(p) (and chosen to fit the detailed 
model described in the next section for n = 6). Notice that the 
activation is concerted and redundant, in that a minimal threshold 
of phosphorylation is required for activation, and activation is 
reached for less than full phosphorylation. 

The MF model reproduces the dynamics of detailed 
allosteric systems 

A validation of the performance of this model for n = 6 is now 
shown in Figure 2 in the context of a system similar to the classical 
and widely used Monod-Wyman-Changeux model of an allosteric 
multisite protein [18]. The original MWC model describes the 
binding of oxygen to the different sites of hemoglobin and the 
allosteric transitions of this protein between two different states. 
Rather than modeling oxygen binding, a protein with i 
phosphorylations is assumed to change between an active 
conformation At and an inactive conformation Each of these 
forms can also be phosphorylated or dephosphorylated at the rates 
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Figure 2. Comparison of MWC and IMF models. A: The MF approximation formula is used to relate the fraction of phosphorylated protein p with 
the active protein concentration P. B: A detailed phosphorylation model structurally similar to the Monod-Wyman-Changeux model [18] is used to 
validate the MF approximation. C: The full MWC model is compared with the MF approximation at steady state as a function of kinase concentration 
E. Also for comparison, a model in which it is simply assumed that the protein has a single site. Here n = 6, £ = 0.1, P = 12 nM, F= 1 nM, L\ =e~"I 2 , 
a. = fi = L.2 = \ for the detailed model, and a= \fi, fl = l, y= — nine, d=L\ for the MF approximation. D: Comparison of the error between the MF and 
the MWC model, for different values of n and e and the remaining parameters computed as above. E: Calculation of the Hill exponent of the MF 
model for different values of n, b. In each case L\ =er"l 2 , to allow for half maximal activation with n/2 phosphorylations, and a=jt=Z2 = 1. F: The 
default equations for the MF approximation. The active protein concentration P is a function of the fraction of active sites p and the total protein P. 
The values of p are calculated via a simple chemical reaction, and a form for the function h(x) is suggested. 
doi:1 0.1 371 /journal.pcbi.1 003443.g002 



given in the diagram in Figure 2B. Although the model is 
interpreted in a different way from MWC, from a mathematical 
point of view it is almost identical. The coefficient e < 1 accounts 
for an assumption that the protein is phosphorylated at a faster 
rate when it is active than when it is inactive. The coefficients 
n,n — 1 etc represent the fact that this is still a nonsequential model: 
for instance, Aq can be phosphorylated at n different sites, so the 
phosphorylation rate is multiplied by n. See the derivation of this 
model from more basic principles in Text S 1 , and a recent review 
on multisite systems by one of the authors [28] . 

Using this multisite model, we now derive parameters for a 
corresponding MF model. For instance, one can define the 'average' 
phosphorylation rate at a given phosphorylation site, regardless of 
whether the protein is active or inactive, !i = ay / E, and the 
dephosphorylation rate /? = /?. By way of derivation of the activity 
function h(x), suppose that a protein with i phosphorylations is 
switching between active and inactive form, 



L x ?< 
Ai ^ /,. 

L 2 

At steady state, we assume that this exchange is balanced and 
calculate L\tfAi = L,%Ii, Then the fraction of active sites with i 
phosphorylations at steady state is 

A Ll B -t e -«(^4 

Ml n) = = : = ; : = r . 

Ii+Ai Lie'+L 2 Li/Li + e-' Li/Lj + e-"^^ 

In other words h(x) = e yx /(5 + e yx ), where )'=— «lns and 
6 = L1/L2. See Text SI. 3 for more details. In particular, the 
ultrasensitive behavior of the function generally increases with the 
number of sites. 
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At any given time, the active protein concentration of the full 
model is defined as Pmwc =^4o+ ■■■+A n . In Figure 2C we 
compare the full 12-variable MWC model for n = 6 with the 
corresponding MF approximation. For every value of the input 
kinase concentration, the resulting concentration of active protein 
P is plotted at steady state. Notice the close similarity between the 
two graphs, which is even more surprising since MF is essentially a 
one-variable model. 

For comparison, we also plot the behavior of an overly 
simplified but all too commonly used model, in which the 
substrate is assumed to have a single phosphorylation site instead 
of n sites, and it is modeled according to the reaction 

Pin <— - ^single site 
HF 

using linear reaction rates. Notice that the behavior of this single- 
site model in Figure 2C is very different from that of the MWC 
model, and that any switch-like behavior in the response is lost. 
This can have important consequences regarding the existence of 
multiple steady states, hysteresis, oscillations etc in the context of 
larger systems, which will be illustrated below. It is easy to show 
that P s i n gie site ~ P P, i.e. it corresponds to Pmf when h(x) = x. It 
should be noted that if the single site system is modeled using 
Michaelis-Menten reactions rather than linear rates, it could have 
strongly ultrasensitive behavior in the saturation regime via zero- 
order ultrasensitivity [33]; see the Discussion section for more 
details. 

We carried out a calculation of the distance between Pmf and 
Pmwc for many different combinations of the number of sites n 
and the allosteric parameter e. For every such set of parameters, 
the two graphs were plotted at steady state as a function of E, and 
the error ma.x e\Pmwc(E) — Pmf(E)\/P was calculated in 
Figure 2D. Notice that the approximation is within 1% precision 
for arbitrary n and s>0.2. On the other hand, in order to obtain 
high ultrasensitivity it is required that n be relatively large and/ or £ 
be small (Figure 2E). See also Figure S2, where additional 
parameter variations are explored over four orders of magnitude 
using the same type of graphs, with similar results. 

It is worth comparing this methodology with the approach 
known in the literature as rule-based modeling, where a series of 
chemical reactions is defined using a streamlined algorithm, and 
high-powered computing is used to handle the resulting large 
number of variables; see e.g. BioNetGen [34]. The advantage of 
this method is that a large number of reactions can be denned and 
handled this way, including complex parameter optimizations. 
One disadvantage is that the combinatorial explosion resulting 
from combining reactions can sometimes exceed the computa- 
tional power. Another is that the large number of equations makes 
any mathematical analysis difficult, if at all possible. 

It is interesting that the MWC model can actually be described 
in terms of rule-based modeling. In Figure SI A we describe a 
series of chemical reactions, over all possible phosphoform states, 
and we show in Text SI that this system is in fact equivalent to the 
MWC model. Thus MF can also be seen as the 1 -variable 
reduction of a system with 2" +1 variables and a much larger 
number of reactions. 

Control theoretic analysis of the pathway reveals two 
robust switches 

In this section, we will embed the MF system within increasingly 
complex systems of equations. We consistendy use upper case for 
proteins and lower case for modified fractions of sites. However, 



we will first provide some technical experimental background 
regarding this specific pathway. 

Background: The cell size checkpoint pathway. Any 

theory for cell size checkpoints should ideally account for several 
mechanistic features. First, cell size checkpoints should translate 
growth into a proportional checkpoint signal. Second, they should 
read the signal to detect when it reaches a threshold that indicates 
when sufficient growth has occurred. Finally, when the threshold is 
reached the checkpoint should trigger a switch-like cell cycle 
transition. All of these mechanisms should be robust and adaptable 
to function in cells of diverse size and shape, and under conditions 
of fluctuating growth rates. 

It was recendy hypothesized that the timing of the polar to 
isotropic growth transition described in the introduction is 
controlled by a mitotic cell size checkpoint. Mitotic cell size 
checkpoints are controlled by Weel and Cdc25 [25]. Weel is a 
protein kinase that delays entry into mitosis by phosphorylating 
and inhibiting Cdkl, while Cdc25 is a phosphatase that promotes 
entry into mitosis by dephosphorylating Cdkl. The budding yeast 
homologs of Weel and Cdc25 are known as Swel and Mihl; 
however, for clarity we will use their more commonly known 
names. In this model, a checkpoint signal originates at the site of 
polar membrane growth in the bud, and downstream components 
read the strength of the signal and trigger activation of Cdkl when 
it reaches a threshold, using the pathway described in Figure IB. A 
switch-like increase in Cdkl during entry into mitosis is postulated 
based on evidence of this qualitative behavior of Cdkl in other 
systems. Also, blocking polar membrane growth causes a Weel- 
dependent arrest before mitosis, which indicates that entry into 
mitosis is linked to membrane growth [25]. The timing of the 
switch from polar to isotropic growth determines the duration of 
polar growth, which influences cell size and shape. In this sense, it 
is a cell size checkpoint. 

Signaling is thought to be initiated by delivery of the GTPase 
Rhol to the site of membrane growth. Rhol is delivered by 
vesicles that are sent to the bud and fuse with the membrane to 
increase its size [35]. Therefore the inflow of Rhol is proportional 
with the rate of polar membrane growth. Rhol on vesicles is 
inactive and it undergoes activation when vesicles fuse at the site of 
growth [35]. Active Rhol binds protein kinase C (Pkcl), a member 
of the protein kinase C family, and induces it to undergo 
autophosphorylation [36]. Importantly, phosphorylation of Pkcl is 
dependent upon and appears to be proportional to membrane 
growth, which suggests that Rhol relays proportional checkpoint 
signals regarding the extent of membrane growth [25]. This 
scenario describes the upstream component of the pathway in 
Figure IB. 

Regarding the downstream component in Figure IB, it is 
thought that hyperphosphorylation of Cdc25 by a poorly 
understood kinase inhibits its activity [37,38]. Thus, a key event 
necessary for initiation of early mitotic events is dephosphorylation 
of Cdc25. This dephosphorylation is carried out by the 
phosphatase PP2A (more precisely, the heterotrimeric complex 
pp2A Cdc55 tha( mcludes the Cdc55 regulatory subunit [25]). PP2A 
is in turn activated by the upstream Pkcl protein. Activation of 
PP2A is also dependent upon an accessory protein named Zdsl 
[25,38] (Zdsl has a redundant paralog named Zds2, which will be 
identified with Zdsl in the model). Zdsl forms a tight stoichio- 
metric complex with PP2A [25]. It also binds to Pkcl, although 
that bond is not included explicitly in the model. PP2A is also 
believed to dephosphorylate Pkcl. Thus, inactivation of PP2A 
causes hyperphosphorylation of Pkcl [25]. Although the functions 
of Pkc 1 phosphorylation are not yet fully known, this observation 
suggests that PP2A restrains activation of Pkcl by checkpoint 
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signals. Similarly, PP2A also dephosphorylates associated Zdsl, 
which suggests that activation of PP2A by Pkcl leads to 
dephosphorylation of Zdsl and activation of targeted dephos- 
phorylation of Cdc25 [38]. In this way PP2A is at the heart of 
checkpoint signaling: it opposes phosphorylation of both Zds 1 and 
Pkcl, and it is responsible for the critical step of dephosphorylating 
Cdc25. There is also good evidence that PP2A regulates Weel, 
although it is less clear whether this function of PP2A is controlled 
by signals from Pkcl [4,38]. 

A bistable switch between Zdsl and PP2A. We can now 
analyze a model involving two different multisite proteins, Zdsl 
and PP2A, each of them described according to the MF 
framework. This forms the downstream component of a larger 
model for the signal transduction of Rhol, and it is described 
inside the lower dashed rectangle in Figure IB. 

The activity of each of the two proteins can be modified through 
phosphorylation and dephosphorylation, and both proteins bind 
together to form a dimer. Only the dimer configuration of these 
molecules is active as an enzyme, and only to the extent that the 
Zdsl and PP2A components have been modified appropriately 
through phosphorylation or dephosphorylation. Moreover, as 
shown in the figure and justified by experimental data, the active 
dimer is itself involved in the modification of the Zdsl protein. 

Denote by Z,P the total monomer concentration of Zdsl and 
PP2A respectively, regardless of their phosphorylation state, Z,P 
their active concentration, z,p the modified fraction of Zdsl and 
PP2A sites, and D the total Zdsl/PP2A dimer concentration. Call 
C the active Rhol /Pkcl dimer concentration, which can be 
treated as a constant input to the system for now. The 
dimerization reaction can be written as 

_ H _ 

Z + P+±D. 

1-2 



r lZ =h 2 ( c W)(i-z), (2) 

where T\ can be calculated from the system parameters. For a 
given value of C this equation can have one or possibly three 
solutions, depending on the parameter values. Plotting the left and 
right hand sides separately can be helpful (Figure 3B). As C 
increases, the function h\(z){\ — z) is rescaled vertically, which can 
control the number of solutions and create a hysteretic response on 
the variable z. 

Given a solution for z in this equation one can solve for the 
concentration of Z, P and every other protein. In particular one 
can calculate the concentration of the active form of the Zds 1 / 
PP2A dimer which we call D, and which constitutes the natural 
output of this system. For specific values of the parameters in the 
model, a bifurcation graph of the output D as a function of the 
input C can be found in Figure 3C (solid line). 

One can interpret this graph in the following way. There is a 
positive feedback loop in the system consisting of Zds 1 promoting 
its own activation via dimerization with PP2A. This positive 
feedback allows the possibility of two stable steady states in the 
system. When C = 0, there can be binding between Zdsl and 
PP2A but since PP2A is inactive, most of the PP2A/Zdsl dimer is 
also inactive. For high values of C, PP2A is forced to become 
active, however the PP2A/Zdsl dimer may still have two steady 
states or only one, depending on the parameters of the system. 

A bistable switch between Rhol and Pkcl. In a similar 
way as it was carried out for the interaction between Zdsl and 
PP2A, one can study the activation of Rho 1 and its binding with 
Pkcl. This upstream system includes all the variables and 
interactions that are not in the Zdsl/PP2A module. Once again 
following the biological evidence and the diagram in Figure 1 B, we 
define the biochemical reactions 



Now, following the argument in the previous section and given 
that both Zdsl and PP2A have been found to have 5 to 10 
phosphorylation sites each, we keep track of z,p using the 
equations 

Z in +± Z, Pin <T± p. 

1-\ 1-3 



16 ?9 V 
R m +±R, (j>+± R in , R -» (j), 

?_6 1-9 ? - 9 

where Ri„,R denote inactive and active Rhol in monomer form 
respectively and v the rate of vesicle fusion to the membrane. That 
is, the protein Rhol is activated at a linear rate, the vesicle flow 
increases Rhol concentration in inactive form, and both active 
and inactive Rhol are degraded. The protein Pkcl is modeled 
using a MF mechanism and denoted by the variables k, K: 



Since Zdsl is activated through dephosphorylation, here z 
represents the fraction of dephosphorylated Zdsl sites, while p 
represents the fraction of phosphorylated PP2A sites. Using the 
MF equation one can then estimate 

Z = Zh i (z), P = Ph 2 (p). 

See Figure 3A for a graph of these two functions. The level of 
ultrasensitivity of each curve is dependent on the estimated 
number of sites in each protein (see the Methods section), but 
notice that the graphs are roughly consistent with the measure- 
ments done for example in the Ste5 protein. A simple mathemat- 
ical analysis of this 5-dimensional model in Text Sl.l shows that 
the solutions of the system converge towards the steady states, and 
that the steady state equations can actually be reduced to a single 
equation for z, 



^ k, K = Kh 3 (k). 

q_ 5 D 

For the study of this subsystem, D will be considered a fixed 
input parameter; it constitutes the feedback from the downstream 
Zdsl/PP2A system. Calling C the total Pkcl/PP2A dimer 
regardless of activity, 

_ 14 _ 

R + K+± C 

1-4 

Notice that Rhol needs to be active to bind with Pkcl. However 
active Rhol binds with Pkcl regardless of the activity or inactivity 
of the latter. Pkc 1 does need to be active in order for the dimer to 
become active as a kinase, so that C / C = K / K = hj,(k), or 
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Figure 3. Control theoretic model analysis. The overall model is analyzed by decomposing it into upstream and downstream levels. A: The 
functions h(x) used to describe the activation of Zds1, PP2A and Pkc1. B: The solutions of the downstream Zdsl - PP2A system (for fixed C) 
correspond to the intersections of the two graphs; see equation (2). C: Bifurcation graph for the downstream Zdsl - PP2A system. D: the solutions of 
the upstream Rho1/Pkc1 system (for fixed D) correspond to the intersections of these graphs; see equation (3). E: Bifurcation graph for the upstream 
Rhol - Pkcl system. F: Both bifurcation graphs superimposed - the steady states of the full model correspond to the intersection of these two graphs. 
doi:1 0.1 371 /journal.pcbi.1 003443.g003 



C=Ch{k). 

Once again, one can carry out a complete mathematical 
analysis of this multi-dimensional system at steady state and 
essentially reduce its solutions to one equation for C. The resulting 
equation is 

where Qi=q-s/qs, 1*2 can be calculated from the parameters 
and constants of the system (other than v, D), and Sk is constant. 
See Text SI. 2 for a full analysis. 

Figure 3D illustrates the solutions of the above equation in the 
general case, plotting right and left hand sides separately. For fixed 
v and D values, C = 0 is always a solution, however there may also 
be solutions C>0. For increasing values of D, the values decrease 
on the right hand side of the equation, and the steady state C > 0 
eventually disappears. Figure 3E illustrates a typical bifurcation 
graph of the steady states of C using the input parameter D and 
frxed v. In the case D = 0 there is actually only one solution for C, 
since equation (3) becomes linear. 

The bifurcation graph can be interpreted as follows. For all but 
small values of the PP2A/Zdsl dimer D, there is the steady state 
with C = k = R = 0, in which the pathway is inactive since the lack 
of active Rhol inhibits the activation of Pkcl. For middle 



concentrations of D, the system may be bistable, again due to the 
positive feedback from Rhol /Pkcl to the activation of Pkcl. For 
high concentrations of D, Pkc 1 is almost fully inactivated by D and 
the pathway is shut down. 

The Rhol pathway can implement a cell size checkpoint 

In order to find the steady states of both subsystems together, 
recall that each one can be reduced to a single equation, so that 
the steady states correspond to the joint solution of the two 
equations. For fixed v, the solutions of the full model form the 
intersection of the graphs for the equations (2), (3). This is 
illustrated in Figure 3F, where the graphs in Figure 3C and 
Figure 3E are superimposed on the same plane. From a control 
perspective, the upstream and downstream systems have each an 
input and an ouptut, and they feed back into each other (see the 
two dotted boxes in Figure IB). The active Zdsl/PP2A dimer D 
also acts as the overall output of the system, since it triggers the 
downstream response to cell cycle regulatory proteins. 

Although it is natural that an increase in the Rhol flow can 
eventually trigger the activation of the pathway, the main focus 
here is not in the flow but in the overall Rho 1 concentration at the 
bud tip. Given that Rhol has a rate of growth proportional to v 
and a linear rate of degradation, at steady state one can show that 
v and R + Rj„ are proportional, R + R in = vqg/q_g. This follows 
from adding the ODE rate equations at steady state, 
0 = R' + R,„' + C = q9V — (Rj„ + R)q_<). In this way one can use 
the total Rhol concentration R = R + Ri„ as a bifurcation 
parameter at steady state even though it is simultaneously a 
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variable in the system. Alternatively, since enzymatic reactions and 
dimer formation are fast processes compared with Rhol flux and 
Rho 1 degradation, one can let R be the slow variable in the system 
and carry out a timescale decomposition analysis using R as a 
constant [39]. 

Let's look at how the system has a hysteretic response for 
increasing values of the flow signal v and the corresponding total 
Rhol concentration at steady state. An increase in these values has 
the effect of raising the graph associated with the upstream system, 
as shown in Figure 4A. For smaller values of R, the intersection of 
both graphs includes three positive steady states (notice the two 
graphs don't quite intersect at the origin). But when R increases 
over a certain threshold, the intersection of the two graphs 
contains a single positive steady state, with a large value of D. This 
can cause an abrupt change in the qualitative behavior of the 
system, triggering a sudden increase in the Zdsl/PP2A output. 
Once this change has taken place, the concentration of the output 
stays high even if the input decreases. 

Figure 4B shows a sample timecourse of the system for a time- 
variable vesicle flow (dotted line). The total Rhol concentration 
increases over time with the inflow of vesicles. At a certain 
timepoint the active Rho 1 concentration abruptly increases, due to 
the switch at the Pkcl/Rhol upstream level. An increase in 



Rhol/Pkcl concentration some time before this can be seen in 
Figure 4C. At a later time the switch between Zdsl and PP2A is 
also triggered, leading to a sudden increase in PP2A/Zdsl 
concentration. Even under variable flow, the total Rhol concen- 
tration roughly corresponds to the membrane accumulated at the 
bud, except for a certain amount of variability due to Rhol 
degradation. Lowering the Rhol degradation rate can decrease 
this difference. Notice that the vesicle flow oscillations do not 
correspond to cell division, but to oscillations in the rate of growth, 
for instance due to varying food availability. 

In Figure 4D we plot the output signal PP2A/Zdsl as a function 
of total Rhol at steady state and overlay the solution of the 
timecourse simulation (red stars). This graph also illustrates the 
hysteretic behavior of the system, in that once a critical threshold 
of Rho 1 concentration is reached, the output signal is dramatically 
increased. This change would constitute a clear signal that the bud 
has reached a large enough size for crossing the polar/isotropic 
growth checkpoint. 

Since both the downstream (PP2A, Zdsl) switch and the 
upstream (Rhol,Pkcl) switch are driven by positive feedback 
loops, it is valid to ask which of the two loops is more relevant for 
the overall system behavior. We argue that it is the downstream 
loop that is more essential, using the bifurcation analysis in Figure 3 
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Figure 4. Checkpoint analysis. A: Qualitative analysis of the steady states of the system for different values of the total Rhol concentration R, as a 
result of changing the value of the fixed flux rate v. When R reaches a sufficiently large value there is a single positive steady state and it has a large D 
concentration. B,C: Response of the checkpoint pathway to a variable vesicle flow input (dashed line). For simplicity the membrane is modeled using 
the equation m'(t) = q 9 -v(t). D: The steady state values of the output D as a function of R (solid line), and plot of the timecourse in 4C over time 
(stars). 

doi:1 0.1 371 /journal.pcbi.1 003443.g004 
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and Figure 4A. If the upstream system is not bistable but has a 
single steady state for every input D, then the graph in Figure 3E is 
replaced by a single-valued decreasing function. Nevertheless this 
(green) line can still have one or three intersections with the (blue) 
downstream dose response in Figure 4A, indicating hysteresis for 
the overall system. On the other hand if the downstream system is 
not bistable, then the blue curve in Figure 4A is replaced by a 
single-valued, increasing function, which would be unlikely to have 
three intersection points with the (green) upstream dose response. 
Thus the downstream switch is essential, while the upstream 
switch is not. 

Notice that this system contains the standard elements of a 
signal transduction pathway, including an initiating signal 
(vesicles), a sensor (Rhol), a series of transducers (Pkcl, PP2A, 
etc), and an effector (active Zds 1 /PP2A). Total Rho 1 is a proxy for 
the membrane concentration, even if bud growth slows for a 
period of time, and the cascade of reactions allows the signal to be 
transduced from the membrane to Zdsl/PP2A and ultimately the 
cyclin dependent kinase. To ensure the high fidelity of the signal 
transmission [40], the downstream signal is sent abruptly after 
total Rhol concentration reaches a particular size. Notice that 
longer periods of inactivity can potentially reduce the Rhol 
concentration significantly - one possible prediction is that after 
such a period the bud grows longer than expected. 

Discussion 

In this paper we have introduced a simple and compact 
framework to describe the dynamics of allosteric multisite 
phosphorylation systems, and we have applied this tool to a new 
molecular model of a size checkpoint in budding yeast. Multisite 
phosphorylation modeling can be problematic because ignoring 
the multiple sites can have significant effects in the dynamics, while 
introducing many auxiliary phosphoform variables can be 
cumbersome in more realistic models. The modified fraction 
approach is intuitive and flexible (model the sites instead of the 
protein), and it only introduces one additional variable per protein. 

The components of the MF model, namely the function h(x) and 
the rates of phosphorylation and dephosphorylation of individual 
sites, can potentially be subject to direct experimental measurement, 
unlike the use of more abstract Hill function terms. This can allow to 
carry out 'raw-data modeling' e.g. to use an experimentally 
measured activity function h(x) directly in the model rather than 
using it to derive parameters. This methodology is also useful out of 
equilibrium when the timescale of phosphorylation is sufficiently fast 
compared with other timescales in the system. 

There are several reasons why the MF method might be 
particularly suitable for modeling many multisite phosphorylation 
systems. Nonsequential phosphorylation is likely more common in 
nature than the more often modeled sequential systems, since 
enforcing sequential phosphorylations would require an additional 
mechanistic effort. Bioinformatic data suggests that most phos- 
phorylation sites in multisite proteins are located in unstructured 
and unconsented protein regions [6], suggesting that often it is the 
collective effect that matters rather than the individual sites. There 
is also experimental evidence in yeast signal transduction that 
certain proteins, such as Ste5, are activated in a concerted and 
redundant manner, although this type of information is still 
unknown for most proteins. Notice that the approximation 
formula would still hold if the protein activation is not concerted 
or redundant. In that case the formula will just approximate a dose 
response that may not be ultrasensitive. 

One of the best known mechanisms for ultrasensitive dose 
responses is zero-order ultrasensitivity, as suggested by Goldbeter 



and Koshland [33,41]. Its main assumption is that substrate 
concentration needs to be in the saturation regime i.e. large 
compared to the K m value of the enzymes. The MF method does 
not pose any constraint on K m , in fact the linear regime we used 
can be found when substrate concentrations are small compared to 
K„, values. Moreover, MF also applies when the enzymatic 
reactions involve complex formation, by writing a Michaelis- 
Menten equation for dpjdt. Therefore zero-order ultrasensitivity 
can be used in synergy with MF in the saturation regime, and MF 
can be used regardless of K m value. A zero-order dose response 
could likely replicate the behavior of the MW C model as shown in 
Figure 2C, however it could not be considered a short hand 
notation for MWC since the two mechanisms are fundamentally 
different. 

In the case of the checkpoint pathway, the active proteins Pkc 1 
and PP2A have been found to have an approximate K m of 0.5 \\.M 
[42] and 1.2 \iM [43], respectively, for specific targets. The 
overall concentrations of their substrates in the cell are much lower 
- however these proteins tend to localize at the bud, so that the 
resulting local concentrations are unknown and it is unclear 
whether a zero-order approach would apply. A recent paper by 
Martins and Swain [44] points out that zero-order ultrasensitivity 
often results from low enzyme to substrate ratios, and localized 
proteins that act as enzymes and substrates for each other would 
likely not satisfy such ratios. That paper proposes instead a 
mechanism involving an allosteric model analogous to MWC, 
using enzyme sequestration to obtain ultrasensitivity. The paper 
by Kapuy et al [32] also proposes a mechanism for bistability 
through ultrasensitive effects, and this mechanism is applied to a 
detailed model of the budding yeast G 1 -S transition in Barik et al 
[45]. Other mechanisms for ultrasensitivity involve competition 
among substrates for the same enzyme [46] and protein 
localization [12], among others [28,47]. 

More generally, in cell regulatory networks there is a need to 
implement nontrivial dynamics such as bistable switches and 
hysteresis, which requires some form of nonlinear response in 
addition to the right feedback interconnections. It has been 
observed that several regulatory proteins have multiple phosphor- 
ylation sites, and there are many open questions regarding their 
intended function. Together with the onerous nature of modeling 
several multisite proteins using sequential networks and multiple 
variables each, it can be seen why a one-variable reduction such as 
MF can allow for much-needed simplicity. 

The actual mechanisms regulating the interactions between cell 
size and cell division remain largely unanswered in many cases. 
This has left few alternative options apart from somewhat heuristic 
approaches in otherwise very detailed models, see e.g. [48]. The 
present model is an attempt, based on recent experiments, to 
construct a detailed mechanistic model in the context of the polar 
to isotropic bud transition in yeast. Notice that if the proteins 
PP2A, Pkcl, Zdsl had only one site each, then hj(x) = x according 
to the argument in the first Results section, and the downstream 
and upstream models could never be bistable (see equations (2) 
and (3)). The multiple sites are providing the underlying 
nonlinearity so that the models can have interesting dynamical 
behaviors. This is consistent with the work by Yang et al [49], 
which reached the same conclusion through randomized param- 
eter searches in multisite cell cycle models. Also, the equations (2), 
(3), which represent the steady states of the downstream and 
upstream systems, have the same qualitative behavior for a wide 
range of parameters. In this sense one can say that the switch-like 
nature of the checkpoint is robust to many parameter changes, 
provided that a few key qualities are satisfied. The bistability in 
each subsystem is due in part to positive feedback loops in each 
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subsystem, one between Pkcl and the Rhol/Pkcl dimer, and 
another between Zdsl and the PP2A/Zdsl dimer. Notice that 
while Rhol/Pkcl activates PP2A, the downstream PP2A/Zdsl 
inactivates Pkcl, forming a negative feedback loop. This feedback 
could serve to reduce the activity of the pathway before a sufficient 
Rhol signal has accumulated. Notice that the switch-like 
activation of Cdkl is a complex process that may well be 
regulated by other mechanisms in conjunction with the switch 
discussed, and that this overall regulation also depends on the 
organism studied. 

The MF framework eliminates several parameters such as the 
number of phosphorylation sites n (as long as it is sufficiendy large), 
the transition rates L, and the cooperativity coefficient e. The 
remaining parameters, such as the shape of the activity function 
h(x), can potentially be measured in the lab using site-directed 
mutagenesis and activity assays. However this is still a formidable 
task and one that is yet to be done for most proteins involved in 
cell cycle regulation. 

Since it is assumed in the derivation of the formula (1) that the 
sites are roughly independent from each other, one might think 
that the MF framework doesn't work for allosteric or cooperative 
systems. However the detailed model in Figure 2 is allosteric, and 
yet the model closely describes its dynamics. In simulations we find 
that the accuracy of the representation is increased when n is large 
(e.g. n> 10) and/or the cooperativity is weak. The use of a MWC- 
type model for multisite phosphorylation has been pointed out in 
the past, see for instance [5] and the more recent [44]. 

Questions for future work include the following: if phosphor- 
ylation and dephosphorylation of the multisite protein P is not 
faster than other processes in the system, can one still approximate 
P(t) away from equilibrium? This might be possible by defining a 
simple differential equation for P instead of the algebraic equation 
(1). Also, the linear dynamics used to calculate the fraction p can 
be replaced by more complex models such as a Michaelis-Menten 
reaction, which may be explored in detail, including the 
interaction with zero-order mechanisms. This might lead to 
bistable behavior in the full multisite model, which raises the 
question of how the corresponding model reduction might be, 
possibly involving multivalued functions h{x). 

Methods 

For convenience we include in one location all chemical 
reactions of the model, mass conservation laws, the definition of 
auxiliary variables following the multisite modeling formalism, and 
a self contained set of differential equations after eliminating 
additional variables. Recall that for multisite proteins x represents 
the fraction of active sites, X the active monomer concentration, 
X the total concentration including active and inactive forms, and 
Sx the total amount of X in the system including dimer and 
monomer forms. Also recall that Rhol and Pkcl are denoted by 
R, K, Zdsl and PP2A by Z, P, and the Rhol/Pkcl and Zdsl/ 
PP2A dimers by C, D, respectively. 

Model reactions 



94 _ ?5 C % 

R+K+±C k in <^ k R in +±R 
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Mass conservation laws 



z,„+z=l Pm+p=\ k in + k=\ 
Z+D = S Z P + D = S P K+C = S K 



MF framework equations 



Z = Zh l (z) = (S z -D)h l (z) 
P=Ph 2 (p) = (Sp-D)h 2 (p) 
K = Kh](k) = (S K -C)h 3 (k) 
C=Ch(k) 
D = bh x (z)h 2 (p) 



Differential equations 
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Parameter values 

Since most quantitative information about the pathway is 
unknown, we make educated estimates on the order of magnitude 
of the parameters. Since parameters are clustered in equations (2) 
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Table 1. Parameter values in the cell size checkpoint model. 
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The model parameters used in the simulation of the checkpoint pathway in Figures 3,4 are detailed in this table, including total protein concentrations, binding and 
unbinding rates, linear reaction rates, rates of vesicle flow and the parameters for the activation functions hj(x). We use the notation Qi = q-i/l<- See the Methods 
section for an estimation of their values. 
doi:1 0.1 371 /journal.pcbi.1 003443.1001 



and (3), dependence on the parameters is more limited. Protein 
concentrations usually range from 0.001 \i.M to 10 (iM in the cell. 
The concentrations of total PP2A (S P ), Zdsl (S z ) and Pkcl (S K ) 
are set between 0.01 (iM and 0.1 \iM as indicated in Table 1. 
Define Q, = q-i/qi for all i. The dissociation rate Q2 has been 
observed to be quite low in experiments since most Zdsl has been 
found bound to PP2A. We set it as 0.001 ixM, which is in the 
range of drugs binding to their target proteins. Q4 is set higher at 
0.1 \iM. The unit-less parameters Q$ , £?6 are set to 0.1 and 1 
respectively, indicating the steady state ratio of inactive to active 
substrate when the two antagonistic enzymes are in similar 
concentration. The rates Qi,Qt, are set to 0.01 [iM, indicating 
that when e.g. C = 0.01 \\M there are equal amounts of active and 
inactive substrate p at steady state. Qg is set at 0.0002 \/s. 

Very little is known about the values of the individual rates 
qi,q-i. Fortunately as it is shown in the analysis in Text SI, most 
of the dynamic rate constants appear only in the form Qi = q-i/qi, 
instead of individually. These steady state ratios are generally 
easier to estimate experimentally than the individual parameters 
[50]. However the actual rates <?i,<7_; determine the transient 
behavior of the system and to some extent determine also its steady 
state values. Since a majority of the reverse rates q-i share the 
same units of 1 / s, we set the values of these parameters and then 
find the corresponding qi to fit the given ratio Qi. For simplicity we 
set q-i =q-2 = <7-3 =<7-4 = (?-6 = 1 S [14,51]. We set 
q_$D= 1 for maximal protein concentration D = Sp, that is, 
q-5 = l/Sp. The Rhol degradation rate q^g is set to 0.0001 s~ l ; 
it can be further decreased in order to stabilize the Rho 1 protein. 

Regarding the activity functions h h we assume that the 
ultrasensitive behavior of these graphs increases with the number 
of phosphorylation sites; see the derivation of h(x) in the Results 
and also [5]. Since PP2A, Zdsl, and Pkcl have been found to have 
around 3, 5, and 8 sites respectively, we implement this with 
parameters that produce the graph observed in Figure 3A. See 
Table 1 for a fist of parameter values. 

The initial conditions used in the model correspond to the 
system in the off state. They are equal to zero for all variables, 
except C= 1 nM and A: = 0.5. 
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Supporting Information 

Figure SI A: A rule-based approach for modeling multisite 
phosphorylation, which results in a system equivalent to MWC. 
Here 7e{0,l}" represents any phosphoform state, and / is the 
result of adding one site phosphorylation to /. B: Timecourse 
simulation for the MW C model in which the kinase concentration 
E is increased linearly over time, and the full MWC model is 
compared with the MF approximation and a model using a single 
site. Here «=12, e = 0.3,P= 12 nM, F=\nM, Li=e,-"l 2 , 
a = P = Li = \ for the detailed model, and c. = ^/e, /i = 1 , 
y= — nine, 8 = L\ for the MF approximation. 
(EPS) 

Figure S2 For the MWC model in Figure 2B, a bootstrap 
parameter analysis of the corresponding MF system. Baseline 
parameter values are a = /? = L2 = l, L\=s~"/ 2 . A: The param- 
eters tt,fi,L\ are varied over four orders of magnitude each, and 
the percentage error with the MF model is calculated. B: A similar 
analysis of the Hill exponent of the MF model. 
(EPS) 

Text SI The supplementary material to this paper has three 
sections. In Section Sl.l and SI. 2 we carry out a mathematical 
analysis of the downstream and upstream subsystems of the 
checkpoint model, respectively. Each analysis results in the 
reduction of the system at steady state to a single equation. In 
Section SI. 3 we derive the MWC model using a system of 
equations involving all 2" protein phosphoforms and a rule-based 
system of reactions. 
(PDF) 
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